{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": [],
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/Xornotor/PPGEEC-CompEvolutiva-Atividades/blob/main/C%C3%B3pia_de_NSGA_II_Otimiza%C3%A7%C3%A3o_de_Microrrede_H%C3%ADbrida.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "# **NSGA-II Otimização de Microrrede Híbrida**\n",
        "\n",
        "**Docente:** Prof. Dr. Edmar Egídio Purcino de Souza\n",
        "\n",
        "**Discentes:** Gabriel Lucas Nascimento Silva, Márcio Barros Costa Júnior, André Paiva Conrado Rodrigues"
      ],
      "metadata": {
        "id": "s6PNZN8nDmml"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# === INSTALAÇÃO E IMPORTS ===\n",
        "!pip install pvlib --quiet\n",
        "!pip install deap --quiet\n",
        "\n",
        "from google.colab import drive\n",
        "import random\n",
        "import numpy as np\n",
        "import pandas as pd\n",
        "import matplotlib.pyplot as plt\n",
        "import pvlib\n",
        "from pvlib.location import Location\n",
        "from pvlib.modelchain import ModelChain\n",
        "from pvlib.pvsystem import PVSystem\n",
        "from pvlib.temperature import TEMPERATURE_MODEL_PARAMETERS\n",
        "from deap import base, creator, tools"
      ],
      "metadata": {
        "id": "Wb4XmJ0Xx3qs"
      },
      "execution_count": 1,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "import warnings\n",
        "\n",
        "# Supressing shapely warnings that occur on import of pvfactors\n",
        "warnings.filterwarnings(action='ignore', module='pvfactors')\n",
        "\n",
        "# Definição de localização do sistema\n",
        "location = Location(latitude=-12.709714, longitude=-38.175534, tz='America/Bahia', altitude=400, name='Camacari')\n",
        "\n",
        "# Importa biblioteca dos parâmetros elétricos dos módulos fotovoltaicos e inversores do GitHub\n",
        "cec_modules = pvlib.pvsystem.retrieve_sam(path='https://raw.githubusercontent.com/NREL/SAM/patch/deploy/libraries/CEC%20Modules.csv')\n",
        "sapm_inverters = pvlib.pvsystem.retrieve_sam(path='https://raw.githubusercontent.com/NREL/SAM/patch/deploy/libraries/CEC%20Inverters.csv')\n",
        "\n",
        "# Definir o modelo de módulo, inversor e ângulo\n",
        "module = cec_modules['JA_Solar_JAM72S30_550_MR']\n",
        "inverter = sapm_inverters['BYD_Auto_Industry_Company_Limited__BEG500KTL_U__480V_']\n",
        "angulo = 14\n",
        "\n",
        "# Usar um modelo de temperatura básico com impacto mínimo (por exemplo, NOCT)\n",
        "temperature_model_parameters = TEMPERATURE_MODEL_PARAMETERS[\"sapm\"][\"open_rack_glass_glass\"]\n",
        "\n",
        "# Criar o array da usina solar, definindo albedo e parâmetros de temperatura com impacto mínimo\n",
        "system = PVSystem(\n",
        "    name=\"Camacari\",\n",
        "    surface_tilt=angulo,\n",
        "    surface_azimuth=0,\n",
        "    module_parameters=module,\n",
        "    inverter_parameters=inverter,\n",
        "    modules_per_string=20,\n",
        "    strings_per_inverter=32,\n",
        "    albedo=0.2,  # Configurando a refletância do solo em 20%\n",
        "    temperature_model_parameters=temperature_model_parameters\n",
        ")\n",
        "\n",
        "# Criar o bloco de parametrização da simulação\n",
        "mc = ModelChain(system, location, aoi_model='no_loss', spectral_model='no_loss')\n",
        "\n",
        "# Ler o arquivo com os dados TMY reais\n",
        "tmy = pd.read_csv('https://raw.githubusercontent.com/222luke/COMPEVOL/refs/heads/main/nrel_camacari_timestamp_updated.csv', sep=';', index_col=0)\n",
        "tmy.index = pd.to_datetime(tmy.index, format=\"%d/%m/%Y %H:%M\")\n",
        "\n",
        "# Simular o modelo\n",
        "mc.run_model(tmy)\n",
        "\n",
        "# Calcular a potência em kW\n",
        "pot = mc.results.ac / 1000\n",
        "\n",
        "# Soma total de energia no período\n",
        "soma_energia = pot.sum()\n",
        "print(f\"Soma total de energia no período: {soma_energia} kWh\")\n"
      ],
      "metadata": {
        "id": "kLQ48yxoznJA",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "b2dd257a-1a4c-47c1-f978-a865b748d350"
      },
      "execution_count": 2,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stderr",
          "text": [
            "/usr/local/lib/python3.11/dist-packages/pvlib/pvsystem.py:2174: UserWarning: Original names contain 505 duplicate(s).\n",
            "  warnings.warn('Original names contain %d duplicate(s).' % n_duplicates)\n",
            "/usr/local/lib/python3.11/dist-packages/pvlib/pvsystem.py:2178: UserWarning: Normalized names contain 505 duplicate(s).\n",
            "  warnings.warn(\n"
          ]
        },
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Soma total de energia no período: 545613.1397531832 kWh\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# Dados da tabela de consumo específico (m³/kWh)\n",
        "dados_consumo_especifico = {\n",
        "    \"Potência (kW)\": [5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75],\n",
        "    \"Motor A\": [3.24, 1.92, 1.48, 1.28, 1.15, 1.04, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan],\n",
        "    \"Motor B\": [3.33, 1.85, 1.41, 1.18, 1.06, 0.96, 0.91, 0.85, 0.82, 0.78, 0.77, 0.75, 0.74, 0.72, 0.70],\n",
        "    \"Motor C\": [np.nan, 1.92, 1.48, 1.25, 1.11, 1.02, 0.95, 0.90, 0.85, 0.83, 0.81, 0.78, np.nan, np.nan, np.nan],\n",
        "    \"Motor D\": [4.27, 2.39, 1.75, 1.42, 1.24, 1.10, 1.01, 0.95, 0.89, 0.85, 0.81, 0.79, 0.77, 0.74, 0.73],\n",
        "    \"Motor E\": [3.71, 2.18, 1.69, 1.42, 1.26, 1.18, 1.10, 1.04, 0.99, 0.96, np.nan, np.nan, np.nan, np.nan, np.nan],\n",
        "}\n",
        "\n",
        "# Criar o DataFrame\n",
        "df_consumo = pd.DataFrame(dados_consumo_especifico)\n",
        "\n",
        "\n",
        "# Dados da tabela de consumo dos motores geradores (em m³/h)\n",
        "dados = {\n",
        "    \"Potência (kW)\": [5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75],\n",
        "    \"Motor A\": [16.19, 19.18, 22.21, 25.52, 28.68, 31.16, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan],\n",
        "    \"Motor B\": [16.67, 18.53, 21.15, 23.66, 26.61, 28.84, 31.70, 34.17, 37.08, 39.03, 42.32, 44.70, 48.13, 50.05, 52.77],\n",
        "    \"Motor C\": [16.87, 19.19, 22.21, 25.02, 27.79, 30.63, 33.19, 35.99, 38.44, 41.64, 44.72, 47.06, np.nan, np.nan, np.nan],\n",
        "    \"Motor D\": [21.36, 23.85, 26.26, 28.42, 31.05, 32.96, 35.23, 37.84, 40.21, 42.57, 44.75, 47.52, 49.76, 51.98, 54.54],\n",
        "    \"Motor E\": [18.54, 21.82, 25.37, 28.41, 31.60, 35.27, 38.41, 41.73, 44.65, 48.13, np.nan, np.nan, np.nan, np.nan, np.nan],\n",
        "}\n",
        "\n",
        "# Criar o DataFrame\n",
        "df = pd.DataFrame(dados)\n",
        "\n",
        "\n",
        "# Dados da Eficiência dos motores geradores (%)\n",
        "dados_eficiencia = {\n",
        "    \"Potência (kW)\": [5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75],\n",
        "    \"Motor A\": [5.34, 9.02, 11.69, 13.56, 15.08, 16.66, 20.91, 22.17, 22.98, 24.26, 24.61, 25.42, 25.58, 26.49, 26.97],\n",
        "    \"Motor B\": [5.68, 10.22, 13.43, 16.01, 17.80, 19.70, 21.97, 22.17, 22.98, 22.76, 23.61, 24.12, 25.08, 25.97, 26.97],\n",
        "    \"Motor C\": [9.87, 9.87, 12.79, 15.14, 17.04, 18.55, 19.97, 21.05, 22.17, 22.74, 23.29, 24.15, 25.10, 25.90, 27.00],\n",
        "    \"Motor D\": [3.58, 6.42, 8.75, 10.78, 12.33, 13.43, 15.21, 16.19, 17.14, 17.99, 18.82, 19.34, 20.17, 20.62, 21.06],\n",
        "    \"Motor E\": [5.01, 8.52, 10.99, 13.08, 14.71, 15.84, 16.94, 17.82, 18.73, 19.31, 19.91, 20.42, 20.98, 21.44, 21.90],\n",
        "}\n",
        "\n",
        "df_eficiencia = pd.DataFrame(dados_eficiencia)\n"
      ],
      "metadata": {
        "id": "dTKjQ8OyD2JB"
      },
      "execution_count": 3,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Parâmetros\n",
        "pot_max_kw = 0     # potência sazonal fixa\n",
        "carga_constante_kw = 120  # carga constante 24h\n",
        "dem_min=pot_max_kw+carga_constante_kw\n",
        "# Índice horário de um ano\n",
        "index = pd.date_range(\"2025-01-01\", periods=8760, freq=\"h\")\n",
        "\n",
        "# Inicializa carga constante 24h\n",
        "carga = pd.Series(carga_constante_kw, index=index)\n",
        "\n",
        "# Horários típicos de irrigação (manhã e tarde)\n",
        "horarios_dia = list(range(5, 10)) + list(range(16, 21))\n",
        "\n",
        "# Mascara booleana para horas de irrigação\n",
        "mask_irriga = carga.index.hour.isin(horarios_dia)\n",
        "\n",
        "# Adiciona carga irrigação nos horários definidos\n",
        "carga.loc[mask_irriga] += pot_max_kw"
      ],
      "metadata": {
        "id": "tlBSz2D-NTsO"
      },
      "execution_count": 4,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#Simulação de microrrede\n",
        "def simular_microrrede(pot, carga, capacidade_bat,\n",
        "                       df_consumo,\n",
        "                       pot_motor_variavel,\n",
        "                       nome_motor_fixo='Motor D',\n",
        "                       pot_motor_fixo=None,\n",
        "                       nome_motor_variavel='Motor B',\n",
        "                       pot_max_bat=120, dod_max=0.8,\n",
        "                       eta_charge=0.9, eta_discharge=0.9,\n",
        "                       demanda_contratada=dem_min, preco_gas=1.80,\n",
        "                       ano_inicial=2025):\n",
        "\n",
        "    # === Tarifas e compensação (Lei 14.300) === A4 VERDE\n",
        "    TE = 0.26248\n",
        "    TUSD_FioB = 0.10142\n",
        "    quota_tusd_fiob = [0.45, 0.60, 0.75, 0.90, 0.9]\n",
        "    anos_ativos = [ano_inicial + i for i in range(len(quota_tusd_fiob))]\n",
        "    valor_liquido_ano = {ano: TE + (1 - quota) * TUSD_FioB for ano, quota in zip(anos_ativos, quota_tusd_fiob)}\n",
        "\n",
        "    tarifas = {\n",
        "        'fora_ponta': 0.3639,\n",
        "        'ponta': 3.1038\n",
        "    }\n",
        "\n",
        "    def obter_cec(nome_motor, potencia_uso):\n",
        "        curva = df_consumo[nome_motor].values\n",
        "        potencias = df_consumo[\"Potência (kW)\"].values\n",
        "        return np.interp(potencia_uso, potencias, curva)\n",
        "\n",
        "    if pot_motor_fixo is None:\n",
        "        potencias_disponiveis = df_consumo[\"Potência (kW)\"].values\n",
        "        pot_motor_fixo = potencias_disponiveis[-1]\n",
        "\n",
        "    # Estado inicial da bateria\n",
        "    soc = 0.5 * capacidade_bat\n",
        "    soc_min = capacidade_bat * (1 - dod_max)\n",
        "    soc_max = capacidade_bat\n",
        "\n",
        "    # Inicializa listas\n",
        "    custo_combustivel_fixo = []\n",
        "    custo_combustivel_var = []\n",
        "    custo_combustivel = []\n",
        "\n",
        "    # Inicialização do dataframe\n",
        "    data_df = pd.DataFrame()\n",
        "    data_df[\"pot\"] = pot\n",
        "    data_df[\"carga\"] = carga\n",
        "    data_df[\"periodo_tarifa\"] = data_df.index.hour.map(lambda x: 'ponta' if 18 <= x <= 20 else 'fora_ponta')\n",
        "    data_df[\"noite\"] = data_df.index.hour.isin(list(range(18, 24)) + list(range(0, 6)))\n",
        "    data_df[\"tarifa\"] = data_df[\"periodo_tarifa\"].map(tarifas)\n",
        "\n",
        "    # Solar atende primeiro\n",
        "    data_df[\"solar_para_carga\"] = data_df[['pot', 'carga']].min(axis=1)\n",
        "    data_df[\"carga_restante\"] = data_df['carga'] - data_df['solar_para_carga']\n",
        "    data_df[\"pot_gerada_restante\"] = data_df['pot'] - data_df['solar_para_carga']\n",
        "\n",
        "    # Motores só operam no período noturno (18h–5h59)\n",
        "    data_df[\"supr_gerador_fixo\"] = 0.0\n",
        "    data_df[\"supr_gerador_var\"] = 0.0\n",
        "    data_df.loc[data_df[\"noite\"] == True, \"supr_gerador_fixo\"] = data_df[\"carga_restante\"].clip(upper=pot_motor_fixo)\n",
        "    data_df.loc[data_df[\"noite\"] == True, \"carga_restante\"] -= data_df.loc[data_df[\"noite\"] == True, \"supr_gerador_fixo\"]\n",
        "\n",
        "    # Motor variável atende restante, respeitando faixa potência\n",
        "    potencias_disponiveis = df_consumo[\"Potência (kW)\"].values\n",
        "    pot_min = potencias_disponiveis.min()\n",
        "    pot_max = potencias_disponiveis.max()\n",
        "    carga_restante_idx = data_df.loc[(data_df[\"noite\"] == True) & (data_df[\"carga_restante\"] > 0)].index\n",
        "    data_df.loc[carga_restante_idx, \"supr_gerador_var\"] = data_df.loc[carga_restante_idx, \"carga_restante\"].clip(pot_min, pot_max)\n",
        "    data_df.loc[carga_restante_idx, \"carga_restante\"] -= data_df.loc[carga_restante_idx, \"supr_gerador_var\"]\n",
        "\n",
        "    # Bateria tenta suprir o que resta\n",
        "    carga_restante_idx = data_df.loc[data_df[\"carga_restante\"] > 0].index\n",
        "    pot_restante_idx = data_df.loc[data_df[\"pot_gerada_restante\"] > 0].index\n",
        "    carga_ou_pot_restante_idx = data_df.loc[(data_df[\"carga_restante\"] > 0) | (data_df[\"pot_gerada_restante\"] > 0)].index\n",
        "    data_df[\"energia_suprida_bat\"] = 0.0\n",
        "    data_df.loc[carga_restante_idx, \"soc\"] = soc\n",
        "\n",
        "    carga_res_arr = data_df.loc[carga_ou_pot_restante_idx, \"carga_restante\"].to_numpy()\n",
        "    pot_res_arr = data_df.loc[carga_ou_pot_restante_idx, \"pot_gerada_restante\"].to_numpy()\n",
        "\n",
        "    carga_real_list = []\n",
        "    supr_bat_list = []\n",
        "    soc_list = []\n",
        "\n",
        "    data_df.loc[data_df.index[0], \"soc\"] = soc\n",
        "\n",
        "    for carga_res, pot_res in zip(carga_res_arr, pot_res_arr):\n",
        "        # Bateria\n",
        "        if carga_res > 0:\n",
        "            energia_extraivel = min((soc - soc_min) * eta_discharge, pot_max_bat)\n",
        "            supr_bat = min(carga_res, energia_extraivel)\n",
        "            supr_bat_list.append(supr_bat)\n",
        "            soc -= supr_bat / eta_discharge\n",
        "\n",
        "        # Armazenar excedente solar\n",
        "        if pot_res > 0:\n",
        "            energia_armazenavel = min(pot_res * eta_charge, pot_max_bat * eta_charge)\n",
        "            espaco_disp = soc_max - soc\n",
        "            carga_real = min(energia_armazenavel, espaco_disp)\n",
        "            soc += carga_real\n",
        "            carga_real_list.append(carga_real)\n",
        "\n",
        "        soc_list.append(soc)\n",
        "\n",
        "    # Fill\n",
        "    data_df.loc[carga_ou_pot_restante_idx, \"soc\"] = soc_list\n",
        "    data_df[\"soc\"] = data_df[\"soc\"].ffill()\n",
        "    data_df.loc[carga_restante_idx, \"energia_suprida_bat\"] = supr_bat_list\n",
        "\n",
        "    # Atualização da carga restante e da potência gerada restante\n",
        "    data_df.loc[carga_restante_idx, \"carga_restante\"] -= data_df.loc[carga_restante_idx, \"energia_suprida_bat\"]\n",
        "    data_df.loc[pot_restante_idx, \"pot_gerada_restante\"] -= np.array(carga_real_list) / eta_charge\n",
        "\n",
        "    # Custo de compra da rede\n",
        "    data_df[\"custo_compra_rede\"] = data_df[\"carga_restante\"] * data_df[\"tarifa\"]\n",
        "\n",
        "    # Créditos por injeção de energia solar\n",
        "    data_df[\"valor_credito\"] = data_df.index.year.map(lambda ano: valor_liquido_ano.get(ano, TE))\n",
        "    data_df[\"credito_energia\"] = data_df[\"pot_gerada_restante\"] * data_df[\"valor_credito\"]\n",
        "\n",
        "    # Consumo específico e custo combustível\n",
        "    data_df[\"cec_fixo\"] = data_df[\"supr_gerador_fixo\"].apply(lambda x: obter_cec(nome_motor_fixo, x))\n",
        "    data_df[\"cec_var\"] = data_df[\"supr_gerador_var\"].apply(lambda x: obter_cec(nome_motor_variavel, x))\n",
        "    data_df[\"custo_comb_m3_fixo\"] = data_df[\"supr_gerador_fixo\"] * data_df[\"cec_fixo\"] * preco_gas\n",
        "    data_df[\"custo_comb_m3_var\"] = data_df[\"supr_gerador_var\"] * data_df[\"cec_var\"] * preco_gas\n",
        "    data_df[\"custo_comb\"] = data_df[\"custo_comb_m3_fixo\"] + data_df[\"custo_comb_m3_var\"]\n",
        "\n",
        "    # Demanda contratada e ultrapassagem\n",
        "    demanda_rede_series = data_df[\"carga_restante\"]\n",
        "    custo_demanda_contratada = 38.61 * demanda_contratada * 12\n",
        "    demanda_max_mensal = demanda_rede_series.resample(\"ME\").max()\n",
        "    excedente_kw_mensal = demanda_max_mensal - demanda_contratada\n",
        "    excedente_kw_mensal[excedente_kw_mensal < 0] = 0\n",
        "    custo_ultrapassagem = (excedente_kw_mensal * 77.22).sum()\n",
        "    custo_total_demanda = custo_demanda_contratada + custo_ultrapassagem\n",
        "\n",
        "    return {\n",
        "        \"soc\": data_df[\"soc\"],\n",
        "        \"energia_suprida_bateria\":  data_df[\"energia_suprida_bat\"],\n",
        "        \"energia_comprada_rede\": data_df[\"carga_restante\"],\n",
        "        \"energia_injetada_solar\": data_df[ \"pot_gerada_restante\"],\n",
        "        \"energia_solar_para_carga\": data_df[\"solar_para_carga\"],\n",
        "        \"energia_suprida_gerador_fixo\": data_df[\"supr_gerador_fixo\"],\n",
        "        \"energia_suprida_gerador_var\": data_df[\"supr_gerador_var\"],\n",
        "        \"custo_combustivel_fixo\": data_df[\"custo_comb_m3_fixo\"],\n",
        "        \"custo_combustivel_var\": data_df[\"custo_comb_m3_var\"],\n",
        "        \"custo_combustivel\": data_df[\"custo_comb\"],\n",
        "        \"custo_compra_rede\": data_df[\"custo_compra_rede\"],\n",
        "        \"credito_geracao_excedente\": data_df[\"credito_energia\"],\n",
        "        \"custo_demanda\": custo_total_demanda,\n",
        "        \"capacidade_utilizada\": capacidade_bat,\n",
        "        \"motor_fixo\": nome_motor_fixo,\n",
        "        \"motor_variavel\": nome_motor_variavel,\n",
        "        \"pot_motor_fixo\": pot_motor_fixo,\n",
        "        \"pot_motor_variavel\": pot_motor_variavel,\n",
        "        \"custo_combustivel_anual\": data_df[\"custo_comb\"].sum()\n",
        "    }\n"
      ],
      "metadata": {
        "id": "Hs_w24VPz4f2"
      },
      "execution_count": 30,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "def lcoe_bateria_obj(x, pot, carga, df_consumo):\n",
        "    capacidade = float(x[0])\n",
        "    pot_var = float(x[1])\n",
        "\n",
        "    if capacidade < 10 or capacidade > 500 or pot_var < 5:\n",
        "        return 1e6\n",
        "\n",
        "    nome_motor_fixo = \"Motor D\"\n",
        "    nome_motor_var = \"Motor B\"\n",
        "\n",
        "    resultado = simular_microrrede(pot, carga, capacidade,\n",
        "                                   df_consumo,\n",
        "                                   pot_motor_variavel=pot_var,\n",
        "                                   nome_motor_fixo=nome_motor_fixo,\n",
        "                                   nome_motor_variavel=nome_motor_var)\n",
        "\n",
        "    energia_bat = resultado[\"energia_suprida_bateria\"].sum()\n",
        "    if energia_bat == 0:\n",
        "        return 1e6\n",
        "\n",
        "    capex = capacidade * 2.8\n",
        "    taxa = 0.08\n",
        "    n = 25\n",
        "    vida_bat = 10\n",
        "\n",
        "    fluxo_reposicoes = sum(capex / ((1 + taxa)**ano) for ano in range(vida_bat, n, vida_bat))\n",
        "    opex_anual = 0.02 * capex\n",
        "    fluxo_opex = sum(opex_anual / ((1 + taxa)**ano) for ano in range(1, n + 1))\n",
        "\n",
        "    ultima_sub = max([ano for ano in range(vida_bat, n + 1, vida_bat) if ano <= n], default=0)\n",
        "    vida_restante = vida_bat - (n - ultima_sub)\n",
        "    salvage = capex * (vida_restante / vida_bat) if vida_restante > 0 else 0\n",
        "    valor_presente_salvage = salvage / ((1 + taxa) ** n)\n",
        "\n",
        "    vpl_total = fluxo_reposicoes + fluxo_opex - valor_presente_salvage\n",
        "    crf = (taxa * (1 + taxa)**n) / ((1 + taxa)**n - 1)\n",
        "    custo_anualizado = (vpl_total + capex) * crf\n",
        "    return custo_anualizado / energia_bat\n",
        "\n",
        "\n",
        "def lcoe_total_obj(x, pot, carga, df_consumo):\n",
        "    capacidade = float(x[0])\n",
        "    pot_var = float(x[1])\n",
        "\n",
        "    # Restrições básicas\n",
        "    if capacidade < 10 or capacidade > 500 or pot_var < 5:\n",
        "        return 1e6\n",
        "\n",
        "    nome_motor_fixo = \"Motor D\"\n",
        "    nome_motor_var = \"Motor B\"\n",
        "\n",
        "    # Simula a microrrede\n",
        "    resultado = simular_microrrede(\n",
        "        pot, carga, capacidade,\n",
        "        df_consumo,\n",
        "        pot_motor_variavel=pot_var,\n",
        "        nome_motor_fixo=nome_motor_fixo,\n",
        "        nome_motor_variavel=nome_motor_var\n",
        "    )\n",
        "\n",
        "    # Energia útil: soma da energia que efetivamente atende à carga\n",
        "    '''\n",
        "    energia_util = (\n",
        "        resultado[\"energia_solar_para_carga\"].sum() +      # solar consumida direto na carga\n",
        "        resultado[\"energia_suprida_bateria\"].sum() +       # bateria\n",
        "        resultado[\"energia_comprada_rede\"].sum() +         # energia da rede\n",
        "        resultado[\"energia_suprida_gerador_fixo\"].sum() +  # motor fixo\n",
        "        resultado[\"energia_suprida_gerador_var\"].sum()     # motor variável\n",
        "    )\n",
        "    '''\n",
        "\n",
        "    energia_util = (\n",
        "        resultado[\"energia_solar_para_carga\"].sum() +      # solar consumida direto na carga\n",
        "        resultado[\"energia_suprida_bateria\"].sum() +       # bateria\n",
        "        resultado[\"energia_injetada_solar\"].sum() +        # solar consumida pela rede\n",
        "        resultado[\"energia_suprida_gerador_fixo\"].sum() +  # motor fixo\n",
        "        resultado[\"energia_suprida_gerador_var\"].sum()     # motor variável\n",
        "    )\n",
        "\n",
        "    if energia_util == 0:\n",
        "        return 1e6\n",
        "\n",
        "    # Parâmetros financeiros e físicos\n",
        "    modules_per_string = 20\n",
        "    strings_per_inverter = 32\n",
        "    custo_modulo = 550\n",
        "    fator_ajuste = 1.85\n",
        "\n",
        "    capex_solar = modules_per_string * strings_per_inverter * custo_modulo * fator_ajuste\n",
        "    capex_bat = capacidade * 2.4\n",
        "    capex_total = capex_solar + capex_bat\n",
        "\n",
        "    taxa = 0.08\n",
        "    n = 25\n",
        "    vida_bat = 10\n",
        "\n",
        "    # VPL Reposições descontado\n",
        "    fluxo_reposicoes = sum(\n",
        "        capex_bat / ((1 + taxa) ** ano) for ano in range(vida_bat, n + 1, vida_bat)\n",
        "    )\n",
        "\n",
        "    # VPL OPEX anual fixo e descontado\n",
        "    opex_anual = 0.02 * capex_bat + 0.02 * capex_solar\n",
        "    fluxo_opex = sum(\n",
        "        opex_anual / ((1 + taxa) ** ano) for ano in range(1, n + 1)\n",
        "    )\n",
        "\n",
        "    # Salvage descontado no último ano\n",
        "    ultima_sub = max([ano for ano in range(vida_bat, n + 1, vida_bat) if ano <= n], default=0)\n",
        "    vida_restante = vida_bat - (n - ultima_sub)\n",
        "    salvage = capex_bat * (vida_restante / vida_bat) if vida_restante > 0 else 0\n",
        "    valor_presente_salvage = salvage / ((1 + taxa) ** n)\n",
        "\n",
        "    # VPL custos variáveis (rede, demanda, combustível), todos descontados\n",
        "    fluxo_rede = sum(\n",
        "        (resultado[\"custo_compra_rede\"].sum() - resultado[\"credito_geracao_excedente\"].sum()) / ((1 + taxa) ** ano)\n",
        "        for ano in range(1, n + 1)\n",
        "    )\n",
        "    fluxo_demanda = sum(\n",
        "        resultado[\"custo_demanda\"] / ((1 + taxa) ** ano) for ano in range(1, n + 1)\n",
        "    )\n",
        "    fluxo_combustivel = sum(\n",
        "        resultado[\"custo_combustivel_anual\"] / ((1 + taxa) ** ano) for ano in range(1, n + 1)\n",
        "    )\n",
        "\n",
        "    # VPL total\n",
        "    vpl_total = capex_total + fluxo_reposicoes + fluxo_opex + fluxo_rede + fluxo_demanda + fluxo_combustivel - valor_presente_salvage\n",
        "\n",
        "    # CRF e custo anualizado\n",
        "    crf = (taxa * (1 + taxa)**n) / ((1 + taxa)**n - 1)\n",
        "    custo_anualizado = vpl_total * crf\n",
        "\n",
        "    # Retorna LCOE\n",
        "    return custo_anualizado / energia_util\n",
        "\n",
        "\n",
        "def fracao_ren_obj(x, pot, carga, df_consumo):\n",
        "    capacidade = float(x[0])\n",
        "    pot_var = float(x[1])\n",
        "\n",
        "    if capacidade < 10 or capacidade > 10000 or pot_var < 5:\n",
        "        return 1e6\n",
        "\n",
        "    nome_motor_fixo = \"Motor D\"\n",
        "    nome_motor_var = \"Motor B\"\n",
        "\n",
        "    resultado = simular_microrrede(pot, carga, capacidade,\n",
        "                                   df_consumo,\n",
        "                                   pot_motor_variavel=pot_var,\n",
        "                                   nome_motor_fixo=nome_motor_fixo,\n",
        "                                   nome_motor_variavel=nome_motor_var)\n",
        "\n",
        "    energia_solar_gerada = pot.sum()\n",
        "    energia_excedente = resultado[\"energia_injetada_solar\"].sum()\n",
        "    energia_solar_util = energia_solar_gerada - energia_excedente\n",
        "\n",
        "    energia_motores = resultado[\"energia_suprida_gerador_fixo\"].sum() + resultado[\"energia_suprida_gerador_var\"].sum()\n",
        "    energia_renovavel_util = energia_solar_util + energia_motores\n",
        "    energia_total = carga.sum()\n",
        "\n",
        "    fracao_ren = energia_renovavel_util / energia_total\n",
        "    return 1 - fracao_ren if fracao_ren >= 0 else 1e6\n",
        "\n",
        "\n"
      ],
      "metadata": {
        "id": "7CJsVeOc6dpM"
      },
      "execution_count": 6,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# %% 3 — Hiperparâmetros\n",
        "POP_SIZE   = 12          # múltiplo de 4\n",
        "N_GEN      = 20\n",
        "CX_PROB    = 0.9\n",
        "MUT_PROB   = 0.2\n",
        "ETA_CROSS  = 15.0\n",
        "ETA_MUT    = 20.0\n",
        "\n",
        "# penalidades iniciais (ficam mais fortes a cada geração)\n",
        "LAMBDA_G0  = 10_000.0     # para desigualdades\n",
        "LAMBDA_H0  = 10_000.0     # para igualdade\n",
        "TOL_H      = 0.01         # tolerância da igualdade\n",
        "\n",
        "# limites das variáveis (capacidade_bateria, pot_motor_variavel)\n",
        "BOUNDS = [(10, 500),    # capacidade_bat (kWh)\n",
        "          (5, 75)]      # pot_motor_variavel (kW)\n",
        "LOWER = [b[0] for b in BOUNDS]\n",
        "UPPER = [b[1] for b in BOUNDS]"
      ],
      "metadata": {
        "id": "DMZY44G97T9S"
      },
      "execution_count": 7,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# %% 5 — Tipos DEAP\n",
        "if \"FitnessMin\" not in dir(creator):\n",
        "    creator.create(\"FitnessMin\", base.Fitness, weights=(-1.0, -1.0))\n",
        "if \"Individual\" not in dir(creator):\n",
        "    creator.create(\"Individual\", list, fitness=creator.FitnessMin)\n",
        "\n",
        "# %% 6 — Toolbox\n",
        "toolbox = base.Toolbox()\n",
        "\n",
        "# Registrando atributos com limites distintos\n",
        "toolbox.register(\"attr_capacidade\", random.uniform, *BOUNDS[0])\n",
        "toolbox.register(\"attr_pot_motor\", random.uniform, *BOUNDS[1])\n",
        "\n",
        "# Cria indivíduo com 2 genes: [capacidade_bat, pot_motor_variavel]\n",
        "toolbox.register(\"individual\", tools.initCycle, creator.Individual,\n",
        "                 (toolbox.attr_capacidade, toolbox.attr_pot_motor), n=1)\n",
        "\n",
        "toolbox.register(\"population\", tools.initRepeat, list, toolbox.individual)\n",
        "\n",
        "\n",
        "# %% 4 — Função Objetivo\n",
        "def evaluate(x):\n",
        "    capacidade = x[0]\n",
        "    pot_motor = x[1]\n",
        "    f1 = lcoe_total_obj(x, pot, carga, df_consumo)\n",
        "    f2 = fracao_ren_obj(x, pot, carga, df_consumo)\n",
        "    return f1, f2\n",
        "\n",
        "toolbox.register(\"evaluate\", evaluate)\n",
        "\n",
        "# Crossover e mutação com limites separados para cada gene\n",
        "def bounded_crossover(ind1, ind2, eta):\n",
        "    for i in range(len(ind1)):\n",
        "        tools.cxSimulatedBinaryBounded(ind1, ind2, low=BOUNDS[i][0], up=BOUNDS[i][1], eta=eta)\n",
        "    return ind1, ind2\n",
        "\n",
        "def bounded_mutation(ind, eta, indpb):\n",
        "    for i in range(len(ind)):\n",
        "        tools.mutPolynomialBounded(ind, low=BOUNDS[i][0], up=BOUNDS[i][1], eta=eta, indpb=indpb)\n",
        "    return ind,\n",
        "\n",
        "toolbox.register(\"mate\", tools.cxSimulatedBinaryBounded,\n",
        "                 low=LOWER, up=UPPER, eta=ETA_CROSS)\n",
        "toolbox.register(\"mutate\", tools.mutPolynomialBounded,\n",
        "                 low=LOWER, up=UPPER, eta=ETA_MUT, indpb=0.5)\n",
        "\n",
        "def sanitize(individual):\n",
        "    for i in range(len(individual)):\n",
        "        val = individual[i]\n",
        "        if isinstance(val, complex):  # evita números complexos\n",
        "            val = val.real\n",
        "        individual[i] = float(min(max(val, LOWER[i]), UPPER[i]))\n",
        "    return individual\n",
        "\n",
        "toolbox.register(\"select\", tools.selNSGA2)"
      ],
      "metadata": {
        "id": "X85_EuRk7YPn"
      },
      "execution_count": 8,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# %% 7 — Loop principal\n",
        "def nsga2_run():\n",
        "    pop = toolbox.population(POP_SIZE)\n",
        "    penalty_stats, feasible_rate = [], []\n",
        "\n",
        "    # Avaliar população inicial\n",
        "    i = 0\n",
        "    print(\"\\rInicializando População...\", end='')\n",
        "    for ind in pop:\n",
        "        i += 1\n",
        "        ind.fitness.values = toolbox.evaluate(ind)\n",
        "    pop = toolbox.select(pop, len(pop))\n",
        "\n",
        "    for gen in range(1, N_GEN + 1):\n",
        "        print(f\"\\rGen {gen}/{N_GEN}...\", end='')\n",
        "        # penalidades crescem ao longo das gerações\n",
        "        # Reprodução\n",
        "        offspring = tools.selTournamentDCD(pop, len(pop))\n",
        "        offspring = list(map(toolbox.clone, offspring))\n",
        "\n",
        "        for c1, c2 in zip(offspring[::2], offspring[1::2]):\n",
        "            if random.random() < CX_PROB:\n",
        "                toolbox.mate(c1, c2)\n",
        "                del c1.fitness.values, c2.fitness.values\n",
        "        for m in offspring:\n",
        "            if random.random() < MUT_PROB:\n",
        "                toolbox.mutate(m)\n",
        "                sanitize(m)\n",
        "                del m.fitness.values\n",
        "\n",
        "        # Avaliar apenas indivíduos novos\n",
        "        for ind in (i for i in offspring if not i.fitness.valid):\n",
        "            ind.fitness.values = toolbox.evaluate(ind)\n",
        "            #print(f\"Obj1: {ind.fitness.values[0]:.4f}, Obj2: {ind.fitness.values[1]:.4f}\")\n",
        "            #print(f\"Indiv: {ind}\")\n",
        "        # Seleção elitista NSGA‑II\n",
        "        pop = toolbox.select(pop + offspring, POP_SIZE)\n",
        "\n",
        "        # Estatísticas\n",
        "\n",
        "        #if gen % 1 == 0 or gen == 1 or gen == N_GEN:\n",
        "            #print(f\"Geração {gen:3d}\")\n",
        "\n",
        "    return pop\n"
      ],
      "metadata": {
        "id": "oh34UzB6Xqpg"
      },
      "execution_count": 9,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# %% 8 — Execução\n",
        "pareto_pop = nsga2_run()"
      ],
      "metadata": {
        "id": "5IS4xyCkYY5L",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "3272f728-29dd-436a-a4b1-373353d37620"
      },
      "execution_count": 31,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Gen 20/20..."
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "pareto_pop"
      ],
      "metadata": {
        "id": "V63JUX92bRPZ",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "3251fcaf-833b-4fb9-9985-6d8b1b20c170"
      },
      "execution_count": 32,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "[[498.3034397586448, 54.97081790988168],\n",
              " [498.3034397586448, 54.97081790988168],\n",
              " [10.109647014719688, 55.47383648676826],\n",
              " [10.109647014719688, 55.47383648676826],\n",
              " [260.86907524862465, 53.14162381990213],\n",
              " [198.74632681865447, 54.97081790988168],\n",
              " [369.57187968126715, 54.718628887884485],\n",
              " [56.035144802617, 54.70420975126952],\n",
              " [97.61748978623174, 54.66549371115015],\n",
              " [97.61748978623174, 54.66549371115015],\n",
              " [450.53172357633423, 54.97081790988168],\n",
              " [376.81631664403164, 55.49231582930811]]"
            ]
          },
          "metadata": {},
          "execution_count": 32
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# %% 9 — Plots\n",
        "# Fronteira aproximada\n",
        "F = np.array([ind.fitness.values for ind in pareto_pop])\n",
        "plt.figure(figsize=(6,5))\n",
        "plt.scatter(F[:,0], F[:,1], s=25, c=\"crimson\")\n",
        "plt.title(\"Fronteira de Pareto aproximada\")\n",
        "plt.xlabel(\"$f_1$\"); plt.ylabel(\"$f_2$\")\n",
        "plt.grid(True); plt.tight_layout(); plt.show()\n"
      ],
      "metadata": {
        "id": "DkkME1-mYcXe",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 507
        },
        "outputId": "a71837d0-f7e3-4590-f5a1-ea9aaf36a4ae"
      },
      "execution_count": 33,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 600x500 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAk0AAAHqCAYAAAAZC3qTAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAARZtJREFUeJzt3X1c1fX9//EnIJcCllfgBROvlpqXiTrtAmskFs1oy6uZIjmbm8yMaYbf0swmuszhzCJb2rY0nV1YPzWTSGpNysJYqdnKNE0DvMhQSUB4//5wnHk8H/STwjkHeNxvN2953p/3+ZzX53Xe5tPP53MOPsYYIwAAAFyQr6cLAAAAqAsITQAAADYQmgAAAGwgNAEAANhAaAIAALCB0AQAAGADoQkAAMAGQhMAAIANhCYAAAAbCE1AHZaTkyMfHx/l5OR4rIaHH35YPj4+Hnt9uN/48eMVHR3t6TKcDB48WIMHD/Z0GajnCE1oUJ577jn5+PhY/nrggQfcWktJSYkefvhhjwYeb3T+exQUFKQf//jHSklJUWFhodvrOXTokB5++GHl5+e7/bUBeJdGni4A8IRHHnlE7du3dxrr3r27W2soKSnRnDlzJOmS/4V8ww036Pvvv1dAQEANVuYdqt6j06dP691339VTTz2ljRs3aseOHQoJCXFbHYcOHdKcOXMUHR2t3r17u+11vdkzzzyjyspKT5cBuB2hCQ3SLbfcopiYGFtzT58+rYCAAPn6et+JWV9fXwUFBV10XklJiVuDRk049z361a9+pWbNmmnRokV69dVXNXr06Eveb2VlpcrKymz1rT6ojffe39+/RvcH1BXe97cA4EFV9witXr1aDz74oNq0aaOQkBAVFxdLktauXau+ffsqODhYzZs311133aWDBw867WP8+PEKDQ3VwYMHlZiYqNDQULVo0ULTpk1TRUWFJGnfvn1q0aKFJGnOnDmOS1EPP/ywYz+7d+/WnXfeqaZNmyooKEgxMTF67bXXLOs99xLf4MGD1b17d+Xl5emGG25QSEiIZs6cKUl69dVXlZCQoNatWyswMFAdO3bU3LlzHXVdzLvvvqt+/fopKChIHTt21NNPP13t3Oeff97Rq6ZNm2rUqFE6cOCArdexctNNN0mS9u7dK0lauHChBg0apGbNmik4OFh9+/bViy++6PI8Hx8fpaSkaOXKlbr66qsVGBioTZs2SZIOHjyou+++WxEREQoMDNTVV1+t5cuXO56bk5Ojfv36SZKSk5Md79Nzzz3nmGNnTVg5duyYpk2bph49eig0NFTh4eG65ZZb9O9//9tpXtV7vGbNGs2cOVORkZFq3Lixhg0b5tLPC733RUVFmjBhgiIiIhQUFKRevXrpr3/9q+O5RUVFatGihQYPHixjjGP8iy++UOPGjTVy5EjH2Pn3NO3bt08+Pj5auHChli5dqg4dOigkJERDhgzRgQMHZIzR3Llz1bZtWwUHB+v222/XsWPHnGr/IWtz2bJl6tixo4KDg9W/f3/985//dJlTVlamWbNmqW/fvmrSpIkaN26s66+/Xlu2bLnoewNUhzNNaJC+++47HTlyxGmsefPmjt/PnTtXAQEBmjZtmkpLSxUQEKDnnntOycnJ6tevn9LT01VYWKjFixfrX//6lz766CNdccUVjudXVFQoPj5eAwYM0MKFC/Xmm2/q8ccfV8eOHfWb3/xGLVq00FNPPaXf/OY3uuOOO/Tzn/9cktSzZ09J0s6dO3XttdeqTZs2euCBB9S4cWP94x//UGJiol566SXdcccdFzy+o0eP6pZbbtGoUaN01113KSIiQtLZ+4VCQ0OVmpqq0NBQvfXWW5o1a5aKi4v12GOPXXCfn3zyiYYMGaIWLVro4Ycf1pkzZzR79mzHvs/1hz/8QQ899JBGjBihX/3qVzp8+LCWLFmiG264waVXdu3Zs0eS1KxZM0nS4sWLNWzYMI0ZM0ZlZWVavXq1hg8frvXr1yshIcHpuW+99Zb+8Y9/KCUlRc2bN1d0dLQKCwv1k5/8xBGqWrRooddff10TJkxQcXGxpk6dqq5du+qRRx7RrFmzdM899+j666+XJA0aNMjRT7tr4nxffvml1q1bp+HDh6t9+/YqLCzU008/rdjYWO3atUutW7d26amPj49mzJihoqIiZWRkKC4uTvn5+QoODnbMs3rvv//+ew0ePFhffPGFUlJS1L59e61du1bjx4/X8ePHde+996ply5Z66qmnNHz4cC1ZskRTpkxRZWWlxo8fr7CwMD355JMXfY9WrlypsrIy/e53v9OxY8f0xz/+USNGjNBNN92knJwczZgxQ1988YWWLFmiadOmOQVUu2vz2Wef1a9//WsNGjRIU6dO1Zdffqlhw4apadOmioqKcswrLi7WX/7yF40ePVoTJ07UiRMn9Oyzzyo+Pl7btm3jUisujQEakBUrVhhJlr+MMWbLli1GkunQoYMpKSlxPK+srMy0bNnSdO/e3Xz//feO8fXr1xtJZtasWY6xpKQkI8k88sgjTq/dp08f07dvX8fjw4cPG0lm9uzZLnX+9Kc/NT169DCnT592jFVWVppBgwaZzp07O8aq6t2yZYtjLDY21kgymZmZLvs995iq/PrXvzYhISFOr2UlMTHRBAUFma+++soxtmvXLuPn52fO/V/Jvn37jJ+fn/nDH/7g9PxPPvnENGrUyGX8fFXv0ZtvvmkOHz5sDhw4YFavXm2aNWtmgoODzddff215LGVlZaZ79+7mpptuchqXZHx9fc3OnTudxidMmGBatWpljhw54jQ+atQo06RJE8f+P/jgAyPJrFixwuX17K4JK6dPnzYVFRVOY3v37jWBgYFOa6fqPW7Tpo0pLi52jP/jH/8wkszixYsdY9W99xkZGUaSef75553qHzhwoAkNDXXa7+jRo01ISIj5z3/+Yx577DEjyaxbt85pf0lJSaZdu3ZOdUsyLVq0MMePH3eMp6WlGUmmV69epry83Ok1AgICnNacnbVZ1fPevXub0tJSx7xly5YZSSY2NtYxdubMGac5xhjz7bffmoiICHP33Xe7vBZgB5fn0CAtXbpUWVlZTr/OlZSU5PSv9w8//FBFRUX67W9/63QvTEJCgrp06aINGza4vMakSZOcHl9//fX68ssvL1rbsWPH9NZbb2nEiBE6ceKEjhw5oiNHjujo0aOKj4/X559/ftHLP4GBgUpOTnYZP/eYqvZ9/fXXq6SkRLt37652fxUVFXrjjTeUmJioH/3oR47xrl27Kj4+3mnuyy+/rMrKSo0YMcJR+5EjRxQZGanOnTvbvjwSFxenFi1aKCoqSqNGjVJoaKheeeUVtWnTxuVYvv32W3333Xe6/vrrtX37dpd9xcbGqlu3bo7Hxhi99NJL+tnPfiZjjFOd8fHx+u677yz3c65LWRPnCgwMdNwnV1FRoaNHjyo0NFRXXXWV5WuPGzdOYWFhjsd33nmnWrVqpY0bN7rs9/z3fuPGjYqMjHS6F8zf319TpkzRyZMn9fbbbzvGn3jiCTVp0kR33nmnHnroIY0dO1a33377BY+lyvDhw9WkSRPH4wEDBkiS7rrrLjVq1MhpvKyszGkd21mbVT2fNGmS04cfxo8f7/S6kuTn5+eYU1lZqWPHjunMmTOKiYm56HsLVIfLc2iQ+vfvf8Ebwc//ZN1XX30lSbrqqqtc5nbp0kXvvvuu01hQUJDjnqUqV155pb799tuL1vbFF1/IGKOHHnpIDz30kOWcoqIiR3iw0qZNG8tP1O3cuVMPPvig3nrrLcd9WlW+++67avd3+PBhff/99+rcubPLtquuusrpL+7PP/9cxhjLuZL9m4iXLl2qH//4x2rUqJEiIiJ01VVXOd2Mv379ej366KPKz89XaWmpY9zqO6POfz8PHz6s48ePa9myZVq2bJnl6xcVFV2wvh+6Js5XWVmpxYsX68knn9TevXud7t2pugR5rvP76ePjo06dOmnfvn1O41bv/VdffaXOnTu7fJiha9euTsciSU2bNtWf//xnDR8+XBEREfrzn/98weM417mBWpIjyJx72ezc8XP/PNhZm1V1nt8Lf39/dejQwaWev/71r3r88ce1e/dulZeXO8bPXw+AXYQmwMK5/+q9FH5+fpf83KqPck+bNs3lLE6VTp06XXAfVvUfP35csbGxCg8P1yOPPKKOHTsqKChI27dv14wZM2rsI+SVlZXy8fHR66+/btmH0NBQW/u5ULD95z//qWHDhumGG27Qk08+qVatWsnf318rVqzQqlWrXOaf34+qY73rrruUlJRk+RpV95fVlnnz5umhhx7S3Xffrblz56pp06by9fXV1KlTL+u9uNy1K0lvvPGGpLOh5uuvv7Z9D1p16766cfPfG85rY20+//zzGj9+vBITEzV9+nS1bNlSfn5+Sk9Pd9wfB/xQhCbAhnbt2kmSPvvsM8enuKp89tlnju0/RHXfol31L2Z/f3/FxcX94P1WJycnR0ePHtXLL7+sG264wTFe9Wm0C2nRooWCg4P1+eefu2z77LPPnB537NhRxhi1b99eP/7xjy+/cAsvvfSSgoKC9MYbbygwMNAxvmLFClvPb9GihcLCwlRRUXHRHlf3Pl3umnjxxRd144036tlnn3UaP378uNOHEqqc33tjjL744gtb4a5du3b6+OOPVVlZ6XS2qeqy17m1btq0SX/5y190//33a+XKlUpKStL777/vdHmtptldm1V1fv755049Ly8v1969e9WrVy/H2IsvvqgOHTro5ZdfdnoPZ8+eXVuHgQaAe5oAG2JiYtSyZUtlZmY6XQp6/fXX9emnn7p8WsuOqu/OOX78uNN4y5YtNXjwYD399NP65ptvXJ53+PDhH/xa0v/+tW/O+Th5WVmZrU9F+fn5KT4+XuvWrdP+/fsd459++qnjrESVn//85/Lz89OcOXOcXqvqtY8ePXpJ9Z9fj4+Pj9MlrX379mndunW2n/+LX/xCL730knbs2OGy/dweN27cWJLr+3S5a8LPz8+lP2vXrq32frW//e1vOnHihOPxiy++qG+++Ua33HLLBV9Hkm699VYVFBRozZo1jrEzZ85oyZIlCg0NVWxsrOMYf/WrX6l///6aN2+e/vKXv2j79u2aN2/eRV/jcthdmzExMWrRooUyMzNVVlbmGH/uuedc3h+rfb7//vvKzc2t6fLRgHCmCbDB399fCxYsUHJysmJjYzV69GjHx8ujo6N13333/eB9BgcHq1u3blqzZo1+/OMfq2nTpurevbu6d++upUuX6rrrrlOPHj00ceJEdejQQYWFhcrNzdXXX3/t8l0+dgwaNEhXXnmlkpKSNGXKFPn4+Ojvf/+7y1/c1ZkzZ442bdqk66+/Xr/97W8df+leffXV+vjjjx3zOnbsqEcffVRpaWnat2+fEhMTFRYWpr179+qVV17RPffco2nTpv3g+s+VkJCgRYsWaejQofrlL3+poqIiLV26VJ06dXKq5ULmz5+vLVu2aMCAAZo4caK6deumY8eOafv27XrzzTcd3yPUsWNHXXHFFcrMzFRYWJgaN26sAQMGqH379pe1Jm677TY98sgjSk5O1qBBg/TJJ59o5cqVlvfmSGfvNbruuuuUnJyswsJCZWRkqFOnTpo4ceJFj/Wee+7R008/rfHjxysvL0/R0dF68cUX9a9//UsZGRmOG8zvvfdeHT16VG+++ab8/Pw0dOhQ/epXv9Kjjz6q22+/3elMTk2yuzb9/f316KOP6te//rVuuukmjRw5Unv37tWKFStc+nbbbbfp5Zdf1h133KGEhATt3btXmZmZ6tatm06ePFkrx4EGwAOf2AM8purj7B988IHl9qqPd69du9Zy+5o1a0yfPn1MYGCgadq0qRkzZozjI/BVkpKSTOPGjV2eO3v2bHP+H7mtW7eavn37moCAAJevH9izZ48ZN26ciYyMNP7+/qZNmzbmtttuMy+++KJLved/5cDVV19tWf+//vUv85Of/MQEBweb1q1bm/vvv9+88cYbLvuozttvv+2ot0OHDiYzM9PyuIwx5qWXXjLXXXedady4sWncuLHp0qWLmTx5svnss88u+BoXe4+qPPvss6Zz584mMDDQdOnSxaxYscKyFklm8uTJlvsoLCw0kydPNlFRUcbf399ERkaan/70p2bZsmVO81599VXTrVs306hRI5evH7CzJqycPn3a/P73vzetWrUywcHB5tprrzW5ubkmNjbW6aPzVe/xCy+8YNLS0kzLli1NcHCwSUhIcPr6B2Mu/N4XFhaa5ORk07x5cxMQEGB69OjhdByvvvqqkWQef/xxp+cVFxebdu3amV69epmysjJjTPVfOfDYY485Pbe6P09W7/EPWZtPPvmkad++vQkMDDQxMTHmnXfecelbZWWlmTdvnmnXrp0JDAw0ffr0MevXr3epHfghfIyx+c9MAIDb5eTk6MYbb9TatWt15513erocoEHjniYAAAAbCE0AAAA2EJoAAABs4J4mAAAAGzjTBAAAYAOhCQAAwAa+3NKGyspKHTp0SGFhYdX+SAUAAOBexhidOHFCrVu3dvmB1LWB0GTDoUOHXH5KNwAA8A4HDhxQ27Zta/11CE02VP2IgQMHDig8PNzD1bhfeXm5Nm/erCFDhsjf39/T5XgN+uKKnlijL67oiTX6Yq26vhQXFysqKsrx93Rt87rQtHTpUj322GMqKChQr169tGTJEvXv399y7s6dOzVr1izl5eXpq6++0p/+9CdNnTrVZd7Bgwc1Y8YMvf766yopKVGnTp20YsUKxcTE2Kqp6pJceHh4gw1NISEhCg8P5w/xOeiLK3pijb64oifW6Iu1i/XFXbfOeNWN4GvWrFFqaqpmz56t7du3q1evXoqPj1dRUZHl/JKSEnXo0EHz589XZGSk5Zxvv/1W1157rfz9/fX6669r165devzxx3XllVfW5qEAAIB6xqvONC1atEgTJ05UcnKyJCkzM1MbNmzQ8uXL9cADD7jM79evn/r16ydJltslacGCBYqKitKKFSscY+3bt6+F6gEAQH3mNaGprKxMeXl5SktLc4z5+voqLi5Oubm5l7zf1157TfHx8Ro+fLjefvtttWnTRr/97W81ceLEap9TWlqq0tJSx+Pi4mJJZ08PlpeXX3ItdVXVMTfEY78Q+uKKnlijL67oiTX6Yq26vri7T14Tmo4cOaKKigpFREQ4jUdERGj37t2XvN8vv/xSTz31lFJTUzVz5kx98MEHmjJligICApSUlGT5nPT0dM2ZM8dlfPPmzQoJCbnkWuq6rKwsT5fgleiLK3pijb64oifW6Iu18/tSUlLi1tf3mtBUWyorKxUTE6N58+ZJkvr06aMdO3YoMzOz2tCUlpam1NRUx+Oqu/OHDBnSYG8Ez8rK0s0338yNieegL67oiTX64oqeWKMv1qrrS9WVIHfxmtDUvHlz+fn5qbCw0Gm8sLCw2pu87WjVqpW6devmNNa1a1e99NJL1T4nMDBQgYGBLuP+/v4NehE39OOvDn1xRU+s0RdX9MQafbF2fl/c3SOv+fRcQECA+vbtq+zsbMdYZWWlsrOzNXDgwEve77XXXqvPPvvMaew///mP2rVrd8n7BAAADY/XnGmSpNTUVCUlJSkmJkb9+/dXRkaGTp065fg03bhx49SmTRulp6dLOnvz+K5duxy/P3jwoPLz8xUaGqpOnTpJku677z4NGjRI8+bN04gRI7Rt2zYtW7ZMy5Yt88xBAgCAOsmrQtPIkSN1+PBhzZo1SwUFBerdu7c2bdrkuDl8//79Tj9b5tChQ+rTp4/j8cKFC7Vw4ULFxsYqJydH0tmvJXjllVeUlpamRx55RO3bt1dGRobGjBnj1mMDAAB1m1eFJklKSUlRSkqK5baqIFQlOjpaxpiL7vO2227TbbfdVhPlAQCABspr7mkCAADwZoQmAAAAG7zu8lxDcWrLNn274FlVFB6VX0QzXTljghrfaP2DiQEAgOcRmjygaEq6Tryw0fH4zNeFKhjxe4WNSVDLDOufoQcAADyLy3NudmrLNqfAdK4TKzeo5O0P3VwRAACwg9DkZt8uePaC24+lP+OmSgAAwA9BaHKzisKjl7UdAAB4BqHJzfwiml3WdgAA4BmEJje7csaEC25vmjbRTZUAAIAfgtDkZo1v7K+wMQmW28LGJCgkNsbNFQEAADv4ygEPaJnxgELviNOx9Gcc39PUNG0igQkAAC9GaPKQkNgYQhIAAHUIl+cAAABsIDQBAADYQGgCAACwgdAEAABgA6EJAADABkITAACADYQmAAAAGwhNAAAANhCaAAAAbCA0AQAA2EBoAgAAsIHQBAAAYAOhCQAAwAZCEwAAgA2EJgAAABsITQAAADYQmgAAAGwgNAEAANhAaAIAALCB0AQAAGADoQkAAMAGQhMAAIANhCYAAAAbCE0AAAA2EJoAAABsIDQBAADYQGgCAACwgdAEAABgA6EJAADABkITAACADYQmAAAAGwhNAAAANhCaAAAAbGjk6QJQ95XtOaATqzbozIECNYqKVNgvExTQMcrTZQEAUKMITbgsxas26PB9f5R8fCRjJB8fHX/iBbXImKHw0bd6ujwAAGqM112eW7p0qaKjoxUUFKQBAwZo27Zt1c7duXOnfvGLXyg6Olo+Pj7KyMi44L7nz58vHx8fTZ06tWaLbqDK9hw4G5gqK6WKCqf/Hp66QOVffu3pEgEAqDFeFZrWrFmj1NRUzZ49W9u3b1evXr0UHx+voqIiy/klJSXq0KGD5s+fr8jIyAvu+4MPPtDTTz+tnj171kbpDdKJVRvOnmGy4iMVr1zv3oIAAKhFXhWaFi1apIkTJyo5OVndunVTZmamQkJCtHz5csv5/fr102OPPaZRo0YpMDCw2v2ePHlSY8aM0TPPPKMrr7yytspvcM4cKDh7Sc6K+e92AADqCa+5p6msrEx5eXlKS0tzjPn6+iouLk65ubmXte/JkycrISFBcXFxevTRRy86v7S0VKWlpY7HxcXFkqTy8nKVl5dfVi11UdUxuxx7u1aqCA48e0nufH6+UrtW9bpf1falAaMn1uiLK3pijb5Yq64v7u6T14SmI0eOqKKiQhEREU7jERER2r179yXvd/Xq1dq+fbs++OAD289JT0/XnDlzXMY3b96skJCQS66lrsvKynIe6B4pLZ584Sdt3Fh7BXkJl76AnlSDvriiJ9boi7Xz+1JSUuLW1/ea0FQbDhw4oHvvvVdZWVkKCgqy/by0tDSlpqY6HhcXFysqKkpDhgxReHh4bZTq1crLy5WVlaWbb75Z/v7+TttOvLhZRx7IOOfTc5KM1Hz+VIXdOcQj9brLhfrSUNETa/TFFT2xRl+sVdeXqitB7uI1oal58+by8/NTYWGh03hhYeFFb/KuTl5enoqKinTNNdc4xioqKvTOO+/oiSeeUGlpqfz8/FyeFxgYaHmPlL+/f4NexFbH33R0gsIG9FLxyvWO72kKH3Ob/Du09VCV7tfQ14UVemKNvriiJ9boi7Xz++LuHnlNaAoICFDfvn2VnZ2txMRESVJlZaWys7OVkpJySfv86U9/qk8++cRpLDk5WV26dNGMGTMsAxN+OP8ObdXsoUmeLgMAgFrlNaFJklJTU5WUlKSYmBj1799fGRkZOnXqlJKTkyVJ48aNU5s2bZSeni7p7M3ju3btcvz+4MGDys/PV2hoqDp16qSwsDB1797d6TUaN26sZs2auYwDAABciFeFppEjR+rw4cOaNWuWCgoK1Lt3b23atMlxc/j+/fvl6/u/b0k4dOiQ+vTp43i8cOFCLVy4ULGxscrJyXF3+QAAoB7zqtAkSSkpKdVejjs/CEVHR8tU9z1B1SBMAQCAS+FVX24JAADgrQhNAAAANhCaAAAAbCA0AQAA2EBoAgAAsIHQBAAAYAOhCQAAwAZCEwAAgA2EJgAAABsITQAAADYQmgAAAGwgNAEAANhAaAIAALCB0AQAAGADoQkAAMAGQhMAAIANhCYAAAAbCE0AAAA2EJoAAABsIDQBAADYQGgCAACwgdAEAABgA6EJAADAhkaeLgBwp1NbtunbBc+qovCo/CKa6coZE9T4xv6eLgsAUAcQmtBgFE1J14kXNjoen/m6UAUjfq+wMQlqmfGABysDANQFXJ5Dg3BqyzanwHSuEys3qOTtD91cEQCgriE0oUH4dsGzF9x+LP0ZN1UCAKirCE1oECoKj17WdgAACE1oEPwiml3WdgAACE1oEK6cMeGC25umTXRTJQCAuorQhAah8Y39FTYmwXJb2JgEhcTGuLkiAEBdw1cOoMFomfGAQu+I07H0Zxzf09Q0bSKBCQBgC6EJDUpIbAwhCQBwSbg8BwAAYAOhCQAAwAZCEwAAgA2EJgAAABsITQAAADYQmgAAAGwgNAEAANhAaAIAALCB0AQAAGADoQkAAMAGQhMAAIANhCYAAAAbCE0AAAA2EJoAAABs8LrQtHTpUkVHRysoKEgDBgzQtm3bqp27c+dO/eIXv1B0dLR8fHyUkZHhMic9PV39+vVTWFiYWrZsqcTERH322We1eAQAAKA+8qrQtGbNGqWmpmr27Nnavn27evXqpfj4eBUVFVnOLykpUYcOHTR//nxFRkZaznn77bc1efJkvffee8rKylJ5ebmGDBmiU6dO1eahAACAeqaRpws416JFizRx4kQlJydLkjIzM7VhwwYtX75cDzzwgMv8fv36qV+/fpJkuV2SNm3a5PT4ueeeU8uWLZWXl6cbbrihho8AAADUV15zpqmsrEx5eXmKi4tzjPn6+iouLk65ubk19jrfffedJKlp06Y1tk8AAFD/ec2ZpiNHjqiiokIRERFO4xEREdq9e3eNvEZlZaWmTp2qa6+9Vt27d692XmlpqUpLSx2Pi4uLJUnl5eUqLy+vkVrqkqpjbojHfiH0xRU9sUZfXNETa/TFWnV9cXefvCY0ucPkyZO1Y8cOvfvuuxecl56erjlz5riMb968WSEhIbVVntfLysrydAleib64oifW6IsremKNvlg7vy8lJSVufX2vCU3NmzeXn5+fCgsLncYLCwurvcn7h0hJSdH69ev1zjvvqG3bthecm5aWptTUVMfj4uJiRUVFaciQIQoPD7/sWuqa8vJyZWVl6eabb5a/v7+ny/Ea9MUVPbFGX1zRE2v0xVp1fam6EuQuXhOaAgIC1LdvX2VnZysxMVHS2ctp2dnZSklJueT9GmP0u9/9Tq+88opycnLUvn37iz4nMDBQgYGBLuP+/v4NehE39OOvDn1xRU+s0RdX9MQafbF2fl/c3SOvCU2SlJqaqqSkJMXExKh///7KyMjQqVOnHJ+mGzdunNq0aaP09HRJZ28e37Vrl+P3Bw8eVH5+vkJDQ9WpUydJZy/JrVq1Sq+++qrCwsJUUFAgSWrSpImCg4M9cJQAAKAu8qrQNHLkSB0+fFizZs1SQUGBevfurU2bNjluDt+/f798ff/3gb9Dhw6pT58+jscLFy7UwoULFRsbq5ycHEnSU089JUkaPHiw02utWLFC48ePr9XjAQAA9YdXhSbp7L1H1V2OqwpCVaKjo2WMueD+LrYdAADADq/5niYAAABvRmgCAACwgdAEAABgA6EJAADABkITAACADYQmAAAAGwhNAAAANhCaAAAAbCA0AQAA2EBoAgAAsIHQBAAAYAOhCQAAwAZCEwAAgA2EJgAAABsITQAAADYQmgAAAGwgNAEAANhAaAIAALCB0AQAAGADoQkAAMAGQhMAAIANhCYAAAAbCE0AAAA2EJoAAABsIDQBAADYQGgCAACwgdAEAABgA6EJAADABkITAACADYQmAAAAGwhNAAAANhCaAAAAbCA0AQAA2EBoAgAAsIHQBAAAYAOhCQAAwAZCEwAAgA2EJgAAABsITQAAADYQmgAAAGwgNAEAANhAaAIAALCB0AQAAGADoQkAAMAGQhMAAIANhCYAAAAbCE0AAAA2EJoAAABs8LrQtHTpUkVHRysoKEgDBgzQtm3bqp27c+dO/eIXv1B0dLR8fHyUkZFx2fsE4F5lew7o6NxMFd7zsI7OzVTZngOeLgkALHlVaFqzZo1SU1M1e/Zsbd++Xb169VJ8fLyKioos55eUlKhDhw6aP3++IiMja2SfANyneNUGHRh0l44vXa2Tr27R8aWrdWDQXSp+YaOnSwMAF14VmhYtWqSJEycqOTlZ3bp1U2ZmpkJCQrR8+XLL+f369dNjjz2mUaNGKTAwsEb2CcA9yvYc0OH7/ihVVkoVFU7/PTx1gcq//NrTJQKAk0aeLqBKWVmZ8vLylJaW5hjz9fVVXFyccnNz3brP0tJSlZaWOh4XFxdLksrLy1VeXn5JtdRlVcfcEI/9QuiLqx/Sk29Xb1RFcODZoHQ+P18de2GDmt5/d02X6BGsFVf0xBp9sVZdX9zdJ68JTUeOHFFFRYUiIiKcxiMiIrR792637jM9PV1z5sxxGd+8ebNCQkIuqZb6ICsry9MleCX64spWT7pHSosnX3jOxvp1mY614oqeWKMv1s7vS0lJiVtf32tCkzdJS0tTamqq43FxcbGioqI0ZMgQhYeHe7AyzygvL1dWVpZuvvlm+fv7e7ocr0FfXP2Qnhz743J998xL1Z5pajLxznp1pom14oyeWKMv1qrrS9WVIHfxmtDUvHlz+fn5qbCw0Gm8sLCw2pu8a2ufgYGBlvdI+fv7N+hF3NCPvzr0xZWdnlw56lad/POqs/cync/XV01HJ9S7vrJWXNETa/TF2vl9cXePvOZG8ICAAPXt21fZ2dmOscrKSmVnZ2vgwIFes08ANSOgY5RaZMyQfH0lP7///tdX8vVVi4wZ8u/Q1tMlAoATrznTJEmpqalKSkpSTEyM+vfvr4yMDJ06dUrJycmSpHHjxqlNmzZKT0+XdPZG7127djl+f/DgQeXn5ys0NFSdOnWytU8AnhM++lYFD+ip4pXrdeZAgRpFRSp8zG0EJgBeyatC08iRI3X48GHNmjVLBQUF6t27tzZt2uS4kXv//v3y9f3fybFDhw6pT58+jscLFy7UwoULFRsbq5ycHFv7BOBZ/h3aqtlDkzxdBgBclFeFJklKSUlRSkqK5baqIFQlOjpaxpjL2icAAIAdXnNPEwAAgDcjNAEAANhAaAIAALCB0AQAAGADoQkAAMAGQhMAAIANhCYAAAAbCE0AAAA2EJoAAABsuOzQ9P333+vgwYMu4zt37rzcXQMAAHiNywpNL774ojp37qyEhAT17NlT77//vmPb2LFjL7s4AAAAb3FZoenRRx9VXl6e8vPztWLFCk2YMEGrVq2SJFs/Ew4AAKCusB2a7r//fp0+fdpprLy8XBEREZKkvn376p133tHTTz+tRx55RD4+PjVbKQAAgAfZDk0ZGRn67rvvJEnjx4/XqVOn1LJlS3388ceOOU2bNlVWVpY+/fRTp3EAAIC6znZoat26tfLz8yVJf//733Xq1Cn9/e9/V8uWLZ3mBQQE6IUXXtDbb79do4UCAAB4ku3Q9Pvf/14/+9nPdP3110uSVq5cqUOHDqlJkyaW86+99tqaqRAAAMAL2A5Nv/vd7/Thhx9q6NChMsZo6dKlGjRokMLDw9W1a1eNGjVK8+fP1+uvv16b9QIAAHhEox8yuWfPnurZs6eee+455ebmqnHjxvr444+Vn5+v/Px8vfrqq/rDH/6gEydO1Fa9AAAAHvGDQlOVzz//3PH7AQMGaMCAAY7HfNUAAACoj2r8x6jwVQMAAKA+4mfPAQAA2EBoAgAAsIHQBAAAYAOhCQAAwAZCEwAAgA2EJgAAABsITQAAADYQmgAAAGwgNAEAANhAaAIAALCB0AQAAGADoQkAAMAGQhMAAIANhCYAAAAbCE0AAAA2EJoAAABsIDQBAADYQGgCAACwoZGnCwCAhqJs70EVr9mkMwcK1CgqUmG/TFBAxyhPlwXAJkITALjJwZsnyq+0XDJG8vHR8SdeUIuMGQoffaunSwNgA5fnAKCWle09ePY3lZVSRYXTfw9PXaDyL7/2bIEAbCE0AUAtO7n2jeo3+kjFK9e7rxgAl4zQBAC17MzBwuo3GunMgQL3FQPgkhGaAKCWNWoTUf1GH6lRVKT7igFwyQhNAFDLQofHV7/RSOFjbnNfMQAuGaEJAGpZQPs2Z3/j6yv5+f33v76Sr69aZMyQf4e2ni0QgC185QAAuEnbN/+i71e/7viepvAxtxGYgDrE6840LV26VNHR0QoKCtKAAQO0bdu2C85fu3atunTpoqCgIPXo0UMbN2502n7y5EmlpKSobdu2Cg4OVrdu3ZSZmVmbhwAAlvyjW6vZQ5MUsexhNXtoEoEJqGO8KjStWbNGqampmj17trZv365evXopPj5eRUVFlvO3bt2q0aNHa8KECfroo4+UmJioxMRE7dixwzEnNTVVmzZt0vPPP69PP/1UU6dOVUpKil577TV3HRYAAKgHvCo0LVq0SBMnTlRycrLjjFBISIiWL19uOX/x4sUaOnSopk+frq5du2ru3Lm65ppr9MQTTzjmbN26VUlJSRo8eLCio6N1zz33qFevXhc9gwUAAHAurwlNZWVlysvLU1xcnGPM19dXcXFxys3NtXxObm6u03xJio+Pd5o/aNAgvfbaazp48KCMMdqyZYv+85//aMiQIbVzIAAAoF7ymhvBjxw5ooqKCkVEOH+fSUREhHbv3m35nIKCAsv5BQX/+6K4JUuW6J577lHbtm3VqFEj+fr66plnntENN9xQbS2lpaUqLS11PC4uLpYklZeXq7y8/AcfW11XdcwN8dgvhL64oifW6IsremKNvlirri/u7pPXhKbasmTJEr333nt67bXX1K5dO73zzjuaPHmyWrdu7XKWqkp6errmzJnjMr5582aFhITUdsleKysry9MleCX64oqeWKMvruiJNfpi7fy+lJSUuPX1vSY0NW/eXH5+fiosdP5xA4WFhYqMtP623MjIyAvO//777zVz5ky98sorSkhIkCT17NlT+fn5WrhwYbWhKS0tTampqY7HxcXFioqK0pAhQxQeHn7Jx1hXlZeXKysrSzfffLP8/f09XY7XoC+u6Ik1+uKKnlijL9aq60vVlSB38ZrQFBAQoL59+yo7O1uJiYmSpMrKSmVnZyslJcXyOQMHDlR2dramTp3qGMvKytLAgQMl/e9ymq+v861bfn5+qqysrLaWwMBABQYGuoz7+/s36EXc0I+/OvTFFT2xRl9c0RNr9MXa+X1xd4+8JjRJZ78eICkpSTExMerfv78yMjJ06tQpJScnS5LGjRunNm3aKD09XZJ07733KjY2Vo8//rgSEhK0evVqffjhh1q2bJkkKTw8XLGxsZo+fbqCg4PVrl07vf322/rb3/6mRYsWeew4AQBA3eNVoWnkyJE6fPiwZs2apYKCAvXu3VubNm1y3Oy9f/9+p7NGgwYN0qpVq/Tggw9q5syZ6ty5s9atW6fu3bs75qxevVppaWkaM2aMjh07pnbt2ukPf/iDJk2a5PbjAwAAdZdXhSZJSklJqfZyXE5OjsvY8OHDNXz48Gr3FxkZqRUrVtRUeQAAoIHymu9pAgAA8GaEJgAAABsITQAAADYQmgAAAGwgNAEAANhAaAIAALCB0AQAAGADoQkAAMAGQhMAAIANhCYAAAAbCE0AAAA2EJoAAABsIDQBAADYQGgCAACwgdAEAABgA6EJAADABkITAACADYQmAAAAGwhNAAAANhCaAAAAbCA0AQAA2EBoAgAAsIHQBAAAYAOhCQAAwAZCEwAAgA2EJgAAABsITQAAADYQmgAAAGwgNAEAANhAaAIAALCB0AQAAGADoQkAAMAGQhMAAIANhCYAAAAbCE0AAAA2EJoAAABsIDQBAADYQGgCAACwgdAEAABgA6EJAADABkITAACADYQmAAAAGwhNAAAANhCaAAAAbCA0AQAA2EBoAgAAsIHQBAAAYAOhCQAAwAavC01Lly5VdHS0goKCNGDAAG3btu2C89euXasuXbooKChIPXr00MaNG13mfPrppxo2bJiaNGmixo0bq1+/ftq/f39tHQIAAKiHvCo0rVmzRqmpqZo9e7a2b9+uXr16KT4+XkVFRZbzt27dqtGjR2vChAn66KOPlJiYqMTERO3YscMxZ8+ePbruuuvUpUsX5eTk6OOPP9ZDDz2koKAgdx0WAACoB7wqNC1atEgTJ05UcnKyunXrpszMTIWEhGj58uWW8xcvXqyhQ4dq+vTp6tq1q+bOnatrrrlGTzzxhGPO//3f/+nWW2/VH//4R/Xp00cdO3bUsGHD1LJlS3cdFgAAqAcaebqAKmVlZcrLy1NaWppjzNfXV3FxccrNzbV8Tm5urlJTU53G4uPjtW7dOklSZWWlNmzYoPvvv1/x8fH66KOP1L59e6WlpSkxMbHaWkpLS1VaWup4XFxcLEkqLy9XeXn5JR5h3VV1zA3x2C+EvriiJ9boiyt6Yo2+WKuuL+7uk9eEpiNHjqiiokIRERFO4xEREdq9e7flcwoKCiznFxQUSJKKiop08uRJzZ8/X48++qgWLFigTZs26ec//7m2bNmi2NhYy/2mp6drzpw5LuObN29WSEjIpRxevZCVleXpErwSfXFFT6zRF1f0xBp9sXZ+X0pKStz6+l4TmmpDZWWlJOn222/XfffdJ0nq3bu3tm7dqszMzGpDU1pamtMZrOLiYkVFRWnIkCEKDw+v/cK9THl5ubKysnTzzTfL39/f0+V4Dfriip5Yoy+u6Ik1+mKtur5UXQlyF68JTc2bN5efn58KCwudxgsLCxUZGWn5nMjIyAvOb968uRo1aqRu3bo5zenatavefffdamsJDAxUYGCgy7i/v3+DXsQN/firQ19c0RNr9MUVPbFGX6yd3xd398hrbgQPCAhQ3759lZ2d7RirrKxUdna2Bg4caPmcgQMHOs2Xzp66q5ofEBCgfv366bPPPnOa85///Eft2rWr4SMAAAD1mdecaZKk1NRUJSUlKSYmRv3791dGRoZOnTql5ORkSdK4cePUpk0bpaenS5LuvfdexcbG6vHHH1dCQoJWr16tDz/8UMuWLXPsc/r06Ro5cqRuuOEG3Xjjjdq0aZP+3//7f8rJyfHEIQIAgDrKq0LTyJEjdfjwYc2aNUsFBQXq3bu3Nm3a5LjZe//+/fL1/d/JsUGDBmnVqlV68MEHNXPmTHXu3Fnr1q1T9+7dHXPuuOMOZWZmKj09XVOmTNFVV12ll156Sdddd53bjw8AANRdXhWaJCklJUUpKSmW26zODg0fPlzDhw+/4D7vvvtu3X333TVRHgAAaKC8LjQBAFBbyvYc0IlVG3TmQIEaRUUq7JcJCugY5emyUEcQmgAADULxqg06fN8fJR8fyRjJx0fHn3hBLTJmKHz0rZ4uD3WA13x6DgCA2lK258DZwFRZKVVUOP338NQFKv/ya0+XiDqA0AQAqPdOrNpw9gyTFR+peOV69xaEOonQBACo984cKDh7Sc6K+e924CK4pwkAUO81ioq84Jkmn7AQHZ2byQ3iuCBCEwCg3gv7ZYKOP/GC9cZKoxPPr5d8fLlBHBfE5TkAQL0X0DFKLTJmSL6+kp/ff//r+7+zT5WGG8RxUZxpAgA0COGjb1XwgJ4qXrnecRmu8vgJFa/ccDYone+/N4g3e2iS+4uFVyI0AQAaDP8ObZ1CUOE9D3ODOGzj8hwAoMG62A3ijaIi3VsQvBqhCQDQYIX9MuGCZ5rCx9zm3oLg1QhNAIAGq9obxH191SJjhvw7tPV0ifAi3NMEAGjQrG4QDx9zG4EJLghNAIAG7/wbxAErXJ4DAACwgdAEAABgA6EJAADABkITAACADYQmAAAAGwhNAAAANhCaAAAAbCA0AQAA2EBoAgAAsIHQBAAAYAOhCQAAwAZCEwAAgA2EJgAAABsITQAAADYQmgAAAGwgNAEAANhAaAIAALCB0AQAAGADoQkAAMAGQhMAAIANhCYAAAAbCE0AAAA2EJoAAABsIDQBAADYQGgCAACwgdAEAABgA6EJAADABkITAACADYQmAAAAGwhNAAAANhCaAAAAbCA0AQAA2OB1oWnp0qWKjo5WUFCQBgwYoG3btl1w/tq1a9WlSxcFBQWpR48e2rhxY7VzJ02aJB8fH2VkZNRw1QAAoL7zqtC0Zs0apaamavbs2dq+fbt69eql+Ph4FRUVWc7funWrRo8erQkTJuijjz5SYmKiEhMTtWPHDpe5r7zyit577z21bt26tg8DAADUQ14VmhYtWqSJEycqOTlZ3bp1U2ZmpkJCQrR8+XLL+YsXL9bQoUM1ffp0de3aVXPnztU111yjJ554wmnewYMH9bvf/U4rV66Uv7+/Ow4FAADUM14TmsrKypSXl6e4uDjHmK+vr+Li4pSbm2v5nNzcXKf5khQfH+80v7KyUmPHjtX06dN19dVX107xAAB4sbI9B3R0bqYK73lYR+dmqmzPAU+XVCc18nQBVY4cOaKKigpFREQ4jUdERGj37t2WzykoKLCcX1BQ4Hi8YMECNWrUSFOmTLFdS2lpqUpLSx2Pi4uLJUnl5eUqLy+3vZ/6ouqYG+KxXwh9cUVPrNEXV/TEWm305cTaN3QkbbHk4yMZI/n46OhfXlLz+VMVdueQGnud2lRdX9y9frwmNNWGvLw8LV68WNu3b5ePj4/t56Wnp2vOnDku45s3b1ZISEhNllinZGVleboEr0RfXNETa/TFFT2xVqN9aSzpzykWG85IF/jwlDc6vy8lJSVufX2vCU3NmzeXn5+fCgsLncYLCwsVGRlp+ZzIyMgLzv/nP/+poqIi/ehHP3Jsr6io0O9//3tlZGRo3759lvtNS0tTamqq43FxcbGioqI0ZMgQhYeHX8rh1Wnl5eXKysrSzTffzD1h56AvruiJNfriip5Yq+m+HPvjcn33zEtSRYXrRj9fNZl4p5ref/dlv05tq64vVVeC3MVrQlNAQID69u2r7OxsJSYmSjp7P1J2drZSUqwSsjRw4EBlZ2dr6tSpjrGsrCwNHDhQkjR27FjLe57Gjh2r5OTkamsJDAxUYGCgy7i/v3+D/sPd0I+/OvTFFT2xRl9c0RNrNdaXr76R3/elUmWl6zZfX+mrb+pU/8/vi7tr95rQJEmpqalKSkpSTEyM+vfvr4yMDJ06dcoRcMaNG6c2bdooPT1dknTvvfcqNjZWjz/+uBISErR69Wp9+OGHWrZsmSSpWbNmatasmdNr+Pv7KzIyUldddZV7Dw4AADdrFBV59l4mKz7/3Q7bvCo0jRw5UocPH9asWbNUUFCg3r17a9OmTY6bvffv3y9f3/994G/QoEFatWqVHnzwQc2cOVOdO3fWunXr1L17d08dAgAAXiPslwk6/sQL1huNFD7mNvcWVMd5VWiSpJSUlGovx+Xk5LiMDR8+XMOHD7e9/+ruYwIAoL4J6BilFhkzdHjqgnM+PSfJSC0yZsi/Q1tPl1ineF1oAgAANSd89K0KHtBTxSvX68yBAjWKilT4mNsITJeA0AQAQD3n36Gtmj00ydNl1Hle843gAAAA3ozQBAAAYAOhCQAAwAZCEwAAgA2EJgAAABsITQAAADYQmgAAAGwgNAEAANhAaAIAALCB0AQAAGADoQkAAMAGQhMAAIANhCYAAAAbCE0AAAA2NPJ0AQAAoH46tWWbvl3wrCoKj8ovopmunDFBjW/s7+myLhmhCQAA1LiiKek68cJGx+MzXxeqYMTvFTYmQS0zHvBgZZeOy3MAAKBGndqyzSkwnevEyg0qeftDN1dUMwhNAACgRn274NkLbj+W/oybKqlZhCYAAFCjKgqPXtZ2b0VoAgAANcovotllbfdWhCYAAFCjrpwx4YLbm6ZNdFMlNYvQBAAAalTjG/srbEyC5bawMQkKiY1xc0U1g68cAAAANa5lxgMKvSNOx9KfcXxPU9O0iXU2MEmEJgAAUEtCYmPqdEg6H5fnAAAAbCA0AQAA2EBoAgAAsIHQBAAAYAOhCQAAwAZCEwAAgA2EJgAAABsITQAAADYQmgAAAGwgNAEAANhAaAIAALCBnz1ngzFGklRcXOzhSjyjvLxcJSUlKi4ulr+/v6fL8Rr0xRU9sUZfXNETa/TFWnV9qfp7uerv6dpGaLLhxIkTkqSoqCgPVwIAAM534sQJNWnSpNZfx8e4K57VYZWVlTp06JDCwsLk4+Pj6XLcrri4WFFRUTpw4IDCw8M9XY7XoC+u6Ik1+uKKnlijL9aq64sxRidOnFDr1q3l61v7dxxxpskGX19ftW3b1tNleFx4eDh/iC3QF1f0xBp9cUVPrNEXa1Z9cccZpircCA4AAGADoQkAAMAGQhMuKjAwULNnz1ZgYKCnS/Eq9MUVPbFGX1zRE2v0xZq39IUbwQEAAGzgTBMAAIANhCYAAAAbCE0AAAA2EJoAAABsIDTVQ0uXLlV0dLSCgoI0YMAAbdu27YLzMzIydNVVVyk4OFhRUVG67777dPr0acf2d955Rz/72c/UunVr+fj4aN26dU7PLy8v14wZM9SjRw81btxYrVu31rhx43To0CGneceOHdOYMWMUHh6uK664QhMmTNDJkydr7LgvxFt7Eh0dLR8fH6df8+fPr7Hjvhh390WSHn74YXXp0kWNGzfWlVdeqbi4OL3//vtOcxrSWpHs9aQhrpVzTZo0ST4+PsrIyHAa9+Rakby3L55cL57oyfjx412Od+jQoU5zamStGNQrq1evNgEBAWb58uVm586dZuLEieaKK64whYWFlvNXrlxpAgMDzcqVK83evXvNG2+8YVq1amXuu+8+x5yNGzea//u//zMvv/yykWReeeUVp30cP37cxMXFmTVr1pjdu3eb3Nxc079/f9O3b1+neUOHDjW9evUy7733nvnnP/9pOnXqZEaPHl3jPTifN/ekXbt25pFHHjHffPON49fJkydrvAdWPNGXqv1kZWWZPXv2mB07dpgJEyaY8PBwU1RU5JjTkNZK1X4u1pOGuFaqvPzyy6ZXr16mdevW5k9/+pPTNk+tFWO8uy+eWi+e6klSUpIZOnSo0/EeO3bMaU5NrBVCUz3Tv39/M3nyZMfjiooK07p1a5Oenm45f/Lkyeamm25yGktNTTXXXnut5fyL/SGusm3bNiPJfPXVV8YYY3bt2mUkmQ8++MAx5/XXXzc+Pj7m4MGDF93f5fDWnhhz9n9s5//Pzl28pS/fffedkWTefPNNYwxrxRjXnhjTcNfK119/bdq0aWN27Njh0gNPrhVjvLcvxnhuvXiqJ0lJSeb222+vtq6aWitcnqtHysrKlJeXp7i4OMeYr6+v4uLilJuba/mcQYMGKS8vz3H69Msvv9TGjRt16623XlYt3333nXx8fHTFFVdIknJzc3XFFVcoJibGMScuLk6+vr4ulyFqkjf3pMr8+fPVrFkz9enTR4899pjOnDlzWa9jh7f0paysTMuWLVOTJk3Uq1cvSawVq55UaWhrpbKyUmPHjtX06dN19dVXu2z31FqRvLsvVdy9Xjz9ZygnJ0ctW7bUVVddpd/85jc6evSoY1tNrRV+YG89cuTIEVVUVCgiIsJpPCIiQrt377Z8zi9/+UsdOXJE1113nYwxOnPmjCZNmqSZM2dech2nT5/WjBkzNHr0aMcPViwoKFDLli2d5jVq1EhNmzZVQUHBJb/WxXhzTyRpypQpuuaaa9S0aVNt3bpVaWlp+uabb7Ro0aJLfi07PN2X9evXa9SoUSopKVGrVq2UlZWl5s2bS2q4a+VCPZEa5lpZsGCBGjVqpClTplhu99Rakby7L5Jn1osnezJ06FD9/Oc/V/v27bVnzx7NnDlTt9xyi3Jzc+Xn51dja4UzTQ1cTk6O5s2bpyeffFLbt2/Xyy+/rA0bNmju3LmXtL/y8nKNGDFCxhg99dRTNVyte7izJ6mpqRo8eLB69uypSZMm6fHHH9eSJUtUWlpaE4dSo2qyLzfeeKPy8/O1detWDR06VCNGjFBRUVEtVF273NmThrZW8vLytHjxYj333HPy8fGpxWrdx519qSvrpab+DI0aNUrDhg1Tjx49lJiYqPXr1+uDDz5QTk5OzRZs+0IevF5paanx8/Nzud47btw4M2zYMMvnXHfddWbatGlOY3//+99NcHCwqaiocJmvC1xjLysrM4mJiaZnz57myJEjTtueffZZc8UVVziNlZeXGz8/P/Pyyy9f5MgunTf3xMqOHTuMJLN79+6Lzr0cnu7L+Tp16mTmzZtnjGm4a+V85/bESn1fK3/605+Mj4+P8fPzc/ySZHx9fU27du2MMZ5bK8Z4d1+suGO9eNufoebNm5vMzExjTM2tFc401SMBAQHq27evsrOzHWOVlZXKzs7WwIEDLZ9TUlIiX1/nZeDn5ydJMj/gxxJWnU35/PPP9eabb6pZs2ZO2wcOHKjjx48rLy/PMfbWW2+psrJSAwYMsP06P5Q398RKfn6+fH19XU4j1zRP9sVKZWWl41/ADXGtWDm3J1bq+1oZO3asPv74Y+Xn5zt+tW7dWtOnT9cbb7whyXNrRfLuvlhxx3rxpj9DX3/9tY4ePapWrVpJqsG1YjteoU5YvXq1CQwMNM8995zZtWuXueeee8wVV1xhCgoKjDHGjB071jzwwAOO+bNnzzZhYWHmhRdeMF9++aXZvHmz6dixoxkxYoRjzokTJ8xHH31kPvroIyPJLFq0yHz00UeOT4GVlZWZYcOGmbZt25r8/Hynj3yWlpY69jN06FDTp08f8/7775t3333XdO7c2W0fI/fGnmzdutX86U9/Mvn5+WbPnj3m+eefNy1atDDjxo2r9Z54qi8nT540aWlpJjc31+zbt898+OGHJjk52QQGBpodO3Y49tOQ1oqdnjTEtWLF6hNhnlorxnhvXzy5XjzRkxMnTphp06aZ3Nxcs3fvXvPmm2+aa665xnTu3NmcPn3asZ+aWCuEpnpoyZIl5kc/+pEJCAgw/fv3N++9955jW2xsrElKSnI8Li8vNw8//LDp2LGjCQoKMlFRUea3v/2t+fbbbx1ztmzZYiS5/Kraz969ey23SzJbtmxx7Ofo0aNm9OjRJjQ01ISHh5vk5GRz4sSJWu7GWd7Yk7y8PDNgwADTpEkTExQUZLp27WrmzZvn9Ie8trm7L99//7254447TOvWrU1AQIBp1aqVGTZsmNm2bZtTXQ1prdjpSUNcK1asQpMn14ox3tkXT68Xd/ekpKTEDBkyxLRo0cL4+/ubdu3amYkTJzqCWpWaWCs+xlzmOWQAAIAGgHuaAAAAbCA0AQAA2EBoAgAAsIHQBAAAYAOhCQAAwAZCEwAAgA2EJgAAABsITQAAADYQmgAAAGwgNAFoEJYsWaJ27dqpUaNGmjZtmqfLAVAH8WNUANR7//73vxUTE6NXX31Vffr0UZMmTRQSEuLpsgDUMY08XQAA1Lb169erf//+uvXWWz1dCoA6jNAEoF7r1KmT9uzZI0ny8fHR2LFj9be//c3DVQGoi7g8B6BeKyoq0sCBA/Wb3/xGd911l0JDQxUaGurpsgDUQdwIDqBeCw0N1b59+3TdddcpMjJSY8eO1ZVXXqk777zT06UBqGMITQDqtY8//liS1KNHD0nSvffey+U5AJeE0ASgXsvPz1enTp3UuHFjSdLgwYMVFhbm4aoA1EWEJgD1Wn5+vnr16uXpMgDUA4QmAPVafn6+evfu7ekyANQDhCYA9VZlZaU++eQTzjQBqBF8TxOAesvX11enTp3ydBkA6gm+pwlAgxIXF6d///vfOnXqlJo2baq1a9dq4MCBni4LQB1AaAIAALCBe5oAAABsIDQBAADYQGgCAACwgdAEAABgA6EJAADABkITAACADYQmAAAAGwhNAAAANhCaAAAAbCA0AQAA2EBoAgAAsIHQBAAAYMP/ByjtqG7TFh2YAAAAAElFTkSuQmCC\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "import pandas as pd\n",
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "from openpyxl import load_workbook\n",
        "from openpyxl.drawing.image import Image as XLImage\n",
        "\n",
        "# === Parâmetros ótimos da solução ===\n",
        "capacidade_otima = 400\n",
        "pot_motor_var = 30\n",
        "\n",
        "# === Parâmetros fixos dos motores ===\n",
        "nome_motor_fixo = \"Motor D\"\n",
        "pot_motor_fixo = 75.0\n",
        "nome_motor_variavel = \"Motor B\"\n",
        "\n",
        "# === Simular microrrede ===\n",
        "resultado_otimo = simular_microrrede(\n",
        "    pot, carga, capacidade_otima,\n",
        "    df_consumo,\n",
        "    pot_motor_variavel=pot_motor_var,\n",
        "    nome_motor_fixo=nome_motor_fixo,\n",
        "    pot_motor_fixo=pot_motor_fixo,\n",
        "    nome_motor_variavel=nome_motor_variavel\n",
        ")\n",
        "\n",
        "# === Parâmetros financeiros ===\n",
        "taxa = 0.08\n",
        "n = 25\n",
        "vida_bat = 10\n",
        "\n",
        "# === CAPEX solar\n",
        "modules_per_string = 20\n",
        "strings_per_inverter = 32\n",
        "capex_solar = modules_per_string * strings_per_inverter * 550 * 1.85\n",
        "capex_bat = capacidade_otima * 2.4\n",
        "capex_total = capex_solar + capex_bat\n",
        "\n",
        "# === Energia útil\n",
        "energia_util = (\n",
        "    resultado_otimo[\"energia_solar_para_carga\"].sum() +\n",
        "    resultado_otimo[\"energia_suprida_bateria\"].sum() +\n",
        "    resultado_otimo[\"energia_comprada_rede\"].sum() +\n",
        "    resultado_otimo[\"energia_suprida_gerador_fixo\"].sum() +\n",
        "    resultado_otimo[\"energia_suprida_gerador_var\"].sum()\n",
        ")\n",
        "\n",
        "# === VPLs\n",
        "vpl_reposicoes = sum(capex_bat / ((1 + taxa) ** ano) for ano in range(vida_bat, n + 1, vida_bat))\n",
        "opex_anual = 0.02 * capex_bat + 0.02 * capex_solar\n",
        "vpl_opex = sum(opex_anual / ((1 + taxa) ** ano) for ano in range(1, n + 1))\n",
        "\n",
        "ultima_sub = max([ano for ano in range(vida_bat, n + 1, vida_bat) if ano <= n], default=0)\n",
        "vida_restante = vida_bat - (n - ultima_sub)\n",
        "salvage = capex_bat * (vida_restante / vida_bat) if vida_restante > 0 else 0\n",
        "valor_presente_salvage = salvage / ((1 + taxa) ** n)\n",
        "\n",
        "vpl_rede = sum(\n",
        "    (resultado_otimo[\"custo_compra_rede\"].sum() - resultado_otimo[\"credito_geracao_excedente\"].sum()) / ((1 + taxa) ** ano)\n",
        "    for ano in range(1, n + 1)\n",
        ")\n",
        "vpl_demanda = sum(resultado_otimo[\"custo_demanda\"] / ((1 + taxa) ** ano) for ano in range(1, n + 1))\n",
        "vpl_combustivel = sum(resultado_otimo[\"custo_combustivel_anual\"] / ((1 + taxa) ** ano) for ano in range(1, n + 1))\n",
        "\n",
        "# === VPL total, LCOE e CRF\n",
        "vpl_total = capex_total + vpl_reposicoes + vpl_opex + vpl_rede + vpl_demanda + vpl_combustivel - valor_presente_salvage\n",
        "crf = (taxa * (1 + taxa) ** n) / ((1 + taxa) ** n - 1)\n",
        "custo_anualizado = vpl_total * crf\n",
        "lcoe_final = custo_anualizado / energia_util\n",
        "\n",
        "# === Fração renovável (conforme imagem)\n",
        "energia_solar_util = resultado_otimo[\"energia_solar_para_carga\"].sum()\n",
        "energia_motores = (\n",
        "    resultado_otimo[\"energia_suprida_gerador_fixo\"].sum() +\n",
        "    resultado_otimo[\"energia_suprida_gerador_var\"].sum()\n",
        ")\n",
        "energia_renovavel_util = energia_solar_util + energia_motores\n",
        "energia_total = carga.sum()\n",
        "fracao_renovavel = energia_renovavel_util / energia_total\n",
        "\n",
        "# === Tabela LCOE + FR\n",
        "df_lcoe = pd.DataFrame({\n",
        "    \"Categoria\": [\n",
        "        \"CAPEX\",\n",
        "        \"VPL Reposições\",\n",
        "        \"VPL OPEX\",\n",
        "        \"VPL Rede (compra - crédito)\",\n",
        "        \"VPL Demanda\",\n",
        "        \"VPL Combustível\",\n",
        "        \"Valor Residual (salvage)\",\n",
        "        \"VPL Total\",\n",
        "        \"Custo Anualizado\",\n",
        "        \"Energia útil total [kWh]\",\n",
        "        \"LCOE [R$/kWh]\",\n",
        "        \"Fração renovável\"\n",
        "    ],\n",
        "    \"Valor [R$]\": [\n",
        "        capex_total,\n",
        "        vpl_reposicoes,\n",
        "        vpl_opex,\n",
        "        vpl_rede,\n",
        "        vpl_demanda,\n",
        "        vpl_combustivel,\n",
        "        -valor_presente_salvage,\n",
        "        vpl_total,\n",
        "        custo_anualizado,\n",
        "        energia_util,\n",
        "        lcoe_final,\n",
        "        fracao_renovavel\n",
        "    ]\n",
        "})\n",
        "\n",
        "# === Fluxo horário detalhado\n",
        "df_fluxos = pd.DataFrame({\n",
        "    \"SOC [%]\": resultado_otimo[\"soc\"]/capacidade_otima,\n",
        "    \"Potência FV [kW]\": pot,\n",
        "    \"Carga [kW]\": carga,\n",
        "    \"Energia FV gerada [kWh]\": pot,\n",
        "    \"Energia solar para carga [kWh]\": resultado_otimo[\"energia_solar_para_carga\"],\n",
        "    \"Energia suprida pela bateria [kWh]\": resultado_otimo[\"energia_suprida_bateria\"],\n",
        "    \"Energia comprada da rede [kWh]\": resultado_otimo[\"energia_comprada_rede\"],\n",
        "    \"Energia injetada da solar [kWh]\": resultado_otimo[\"energia_injetada_solar\"],\n",
        "    \"Energia motor fixo [kWh]\": resultado_otimo[\"energia_suprida_gerador_fixo\"],\n",
        "    \"Energia motor variável [kWh]\": resultado_otimo[\"energia_suprida_gerador_var\"],\n",
        "    \"Custo da compra de energia [R$]\": resultado_otimo[\"custo_compra_rede\"],\n",
        "    \"Custo do combustível motor fixo [R$]\": resultado_otimo[\"custo_combustivel_fixo\"],\n",
        "    \"Custo do combustível motor variável [R$]\": resultado_otimo[\"custo_combustivel_var\"],\n",
        "    \"Crédito por injeção excedente [R$]\": resultado_otimo[\"credito_geracao_excedente\"],\n",
        "}, index=pot.index)\n",
        "\n",
        "df_fluxos[\"Capacidade da bateria [kWh]\"] = \"\"\n",
        "df_fluxos.loc[df_fluxos.index[0], \"Capacidade da bateria [kWh]\"] = capacidade_otima\n",
        "\n",
        "# === Resumo de demanda\n",
        "custo_total_demanda = resultado_otimo[\"custo_demanda\"]\n",
        "custo_demanda_contratada = 38.61 * 120 * 12\n",
        "custo_ultrapassagem = custo_total_demanda - custo_demanda_contratada\n",
        "\n",
        "df_resumo = pd.DataFrame({\n",
        "    \"Custo de demanda contratada [R$/ano]\": [custo_demanda_contratada],\n",
        "    \"Custo de ultrapassagem [R$/ano]\": [custo_ultrapassagem],\n",
        "    \"Custo total de demanda [R$/ano]\": [custo_total_demanda],\n",
        "})\n",
        "\n",
        "# === Gráfico de pizza das fontes\n",
        "energia_fontes = {\n",
        "    \"FV direta\": energia_solar_util,\n",
        "    \"Bateria\": resultado_otimo[\"energia_suprida_bateria\"].sum(),\n",
        "    \"Rede\": resultado_otimo[\"energia_comprada_rede\"].sum(),\n",
        "    \"Gerador fixo\": resultado_otimo[\"energia_suprida_gerador_fixo\"].sum(),\n",
        "    \"Gerador variável\": resultado_otimo[\"energia_suprida_gerador_var\"].sum()\n",
        "}\n",
        "\n",
        "fig, ax = plt.subplots()\n",
        "ax.pie(energia_fontes.values(), labels=energia_fontes.keys(), autopct=\"%1.1f%%\", startangle=90)\n",
        "ax.axis(\"equal\")\n",
        "plt.title(\"Contribuição de Energia por Fonte\")\n",
        "plt.savefig(\"grafico_contribuicao_fontes.png\", dpi=300)\n",
        "plt.close()\n",
        "\n",
        "# === Salvar Excel com imagem\n",
        "with pd.ExcelWriter(\"fluxos_microrrede.xlsx\", engine=\"openpyxl\") as writer:\n",
        "    df_fluxos.to_excel(writer, sheet_name=\"Fluxo Horário\")\n",
        "    df_resumo.to_excel(writer, sheet_name=\"Totais e Custos\", index=False)\n",
        "    df_lcoe.to_excel(writer, sheet_name=\"LCOE Detalhado\", index=False)\n",
        "\n",
        "wb = load_workbook(\"fluxos_microrrede.xlsx\")\n",
        "ws = wb.create_sheet(\"Gráfico de Energia\")\n",
        "img = XLImage(\"grafico_contribuicao_fontes.png\")\n",
        "img.anchor = \"A1\"\n",
        "ws.add_image(img)\n",
        "wb.save(\"fluxos_microrrede.xlsx\")\n",
        "\n",
        "print(\"Planilha 'fluxos_microrrede.xlsx' gerada com gráfico.\")\n",
        "print(f\"Fração renovável: {fracao_renovavel:.2%}\")\n"
      ],
      "metadata": {
        "id": "n2ceqhmZxrA8",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "82df7a8e-458e-4970-be1f-2ef64c957e2b"
      },
      "execution_count": 34,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Planilha 'fluxos_microrrede.xlsx' gerada com gráfico.\n",
            "Fração renovável: 84.06%\n"
          ]
        }
      ]
    }
  ]
}